#######################################################################################
########### Not All Intrastate Violence is Civil War: Varieties of Violence ###########
###########       and the Perils of Model Dependence in Civil War           ###########
#######################################################################################
######################### Noel Anderson & Alec Worsnop ################################
#######################################################################################
############################### 26 April 2012 #########################################
#######################################################################################

#########################################################
#########################################################
############# Libraries, Data, and Packages #############
#########################################################
#########################################################

## Load libraries
library(foreign)
library(MatchIt)
library(cem)
library(Zelig)
library(WhatIf)
library(Amelia)
library(xtable)
library(mvtnorm)
library(MASS)
library(lmtest)
library(lattice)
library(apsrtable)
library(sandwich)
library(lmtest)
library(lme4)
library(fields)
library(car)

## Load data
#repdata <- read.dta("C:\\Documents and Settings\\Noel Anderson\\My Documents\\Dropbox\\[GOV2001] Advanced Quant Methods\\Replication Paper\\Lujala Replication data - data Deadly Combat over Natural Resources JCR 2009.dta")
#repdata <- read.dta("/Users/Alec/Documents/Dropbox/mit/Second Year/Second Semester/Gov 2001/paper/Lujalareplication.dta")

## Function to cluster standard errors
cl <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) 
}

## Function to get robust variance covariance matrix
clv <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  return(dfc*sandwich(fm, meat=crossprod(uj)/N))
  }

######################################################################
######################################################################
############# Section 2 - The Perils of Model Dependence #############
######################################################################
######################################################################

#############################################################
#### Section 2.1 - Counterfactual Analysis - Convex Hull ####
#############################################################

#### M2 (total deaths) ####

## Create Original Dataset
keep.vars <- c("ccode","BTTLDtotLN","DRUGP","ALLGEMP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")
all.data = repdata[keep.vars]
data = na.omit(all.data)

### Create counterfactual dataset
DRUGPCF = 1-repdata$DRUGP
repdataCF = cbind(repdata, DRUGPCF)
keep.vars <- c("ccode","BTTLDtotLN","DRUGPCF","ALLGEMP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")
all.data = repdataCF[keep.vars]
cfdata = na.omit(all.data)

### Perform Test
my.result <- whatif(data=data, cfact=cfdata)
summary(my.result)

### Examine Results 
my.result$in.hull
my.result$geom.var
sort(my.result$sum.stat) #max = 0.192622951; cf#110; Iran(1967), min = 0.004098361; cf#161; China (1946))
median(my.result$sum.stat) #median = 0.02868852; cf#7; Burundi(1965))
mean(my.result$sum.stat)

#### M6 (combat intensity) ####

### Create Original Dataset
keep.vars <- c("ccode","BTTLDintLN","DRUGP","ALLGEMP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM","Duration10")
all.data <- repdata[keep.vars]
data <- na.omit(all.data)

### Create counterfactual dataset (which now also includes 'Duration10'')
DRUGPCF = 1-repdata$DRUGP
repdataCF = cbind(repdata, DRUGPCF)
keep.vars <- c("ccode","BTTLDintLN","DRUGPCF","ALLGEMP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM","Duration10")
all.data = repdataCF[keep.vars]
cfdata = na.omit(all.data)

### Perform Test
my.result2 <- whatif(data=data, cfact=cfdata)
summary(my.result2)

### Examine Results
my.result2$in.hull
my.result2$geom.var #0.1584325
sort(my.result2$sum.stat) #max = 0.225409836; cf#110; Iran(1967), min = 0.004098361; cf#161; China (1946))
median(my.result2$sum.stat) #median = 0.02868852; cf#25; Morocco(1971))
mean(my.result2$sum.stat)

## Plot Cumulative Frequencies of Distances to CFs (Min, Max, Median)
par(mfrow=c(1,2))
plot(my.result, type="f", numcf=(c(161,7,110)))
abline(v=0.1610866, col="red")
text(x=0.5, y=0.18, "China (1967)",cex=0.7)
text(x=0.52, y=0.55, "Burundi (1965)",cex=0.7)
text(x=0.22, y=0.85, "Iran (1967)",cex=0.7)
plot(my.result2, type="f", numcf=(c(161,25,110)))
abline(v=0.1584325, col="red")
text(x=0.48, y=0.18, "China (1967)",cex=0.7)
text(x=0.55, y=0.56, "Morocco (1971)",cex=0.7)
text(x=0.215, y=0.85, "Iran (1967)",cex=0.7)

#############################################################
############# Section 2.2 - Matching Analysis  ##############
#############################################################

##################################
##### Model 2 - Total Deaths #####
##################################

keep.vars <- c("ccode","BTTLDtotLN","DRUGP","ALLGEMP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")
all.data <- repdata[keep.vars]
data <- na.omit(all.data)

coefs=c()
coefsne=c()
tstats=c()
n=c()
pre=matrix(,1000,14)
post=matrix(,1000,14)
meanstat=matrix(,1000,13)
L1s=matrix(,1000,13)
minb=matrix(,1000,13)
maxb=matrix(,1000,13)
tfb=matrix(,1000,13)
ftb=matrix(,1000,13)
sfb=matrix(,1000,13)
maxb=matrix(,1000,13)

for (i in 1:1000){
  genetic <- matchit(DRUGP~ALLGEMP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+DEM, data=data, method="genetic")
  genetic.data <- match.data(genetic)
  
  pre.balance.g <- summary(genetic)$sum.all
  post.balance.g <- summary(genetic)$sum.matched
  
  genetic.imbalance <- imbalance(group=genetic.data$DRUGP, data=genetic.data, drop=c("DRUGP","BTTLDtotLN","ccode","distance","weights"))
  
  pre[i,]=(pre.balance.g[15:2,1]-pre.balance.g[15:2,2])/pre.balance.g[15:2,3]
  post[i,]=(post.balance.g[15:2,1]-post.balance.g[15:2,2])/post.balance.g[15:2,3]
  meanstat[i,]=genetic.imbalance$tab[,1]
  L1s[i,]=genetic.imbalance$tab[,3]
  minb[i,]=genetic.imbalance$tab[,4]
  tfb[i,]=genetic.imbalance$tab[,5]
  ftb[i,]=genetic.imbalance$tab[,6]
  sfb[i,]=genetic.imbalance$tab[,7]
  maxb[i,]=genetic.imbalance$tab[,8]
  
  M2 <- lm(BTTLDtotLN~ALLGEMP+DRUGP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+DEM, data=genetic.data)
  M2$se=abs(round(cl(genetic.data, M2, genetic.data$ccode), 2))[,3]
  coefsne=c(coefsne,M2$coefficients[3])
  coefs=c(coefs,round(exp(M2$coefficients),3)[3])
  tstats=c(tstats,M2$se[3])
  n=c(n,dim(genetic.data)[1])
  print(i)
}

## Code for plot (see below)
pre[is.infinite(pre)] <- NA 
post[is.infinite(post)] <- NA

## Code for table in APPENDIX A
pre.imbalance <- imbalance(group=data$DRUGP, data=data, drop=c("DRUGP","BTTLDtotLN","ccode","distance","weights"))
apply(maxb,2,median)
balance2=cbind(pre.imbalance$tab[,1],apply(meanstat,2,mean),pre.imbalance$tab[,3],apply(L1s,2,mean),
pre.imbalance$tab[,5],apply(tfb,2,mean),pre.imbalance$tab[,6],apply(ftb,2,mean),pre.imbalance$tab[,7],apply(sfb,2,mean))
rownames(balance2) = c("Gem Production (In Zone)","Hydrocarbon Prod. (In Zone)","Hydrocarbon Prod. (Out Zone)","Mountainous Terrain","Conflict at Border","Weak Rebel Group","Equal Strength","Cold War", "Internationalized Conflict","Population","Ethnic Fractionalization","Ethnic Polarization","Democracy")
print(xtable(balance2, include.rownames=F,digit=4))

######################################
##### Model 6 - Combat Intensity #####
######################################

keep.vars <- c("ccode","BTTLDintLN","DRUGP","ALLGEMP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM","Duration10")
all.data <- repdata[keep.vars]
data <- na.omit(all.data)

coefs2=c()
coefsne2=c()
tstats2=c()
n2=c()
pre2=matrix(,1000,14)
post2=matrix(,1000,14)
meanstat2=matrix(,1000,14)
L1s2=matrix(,1000,14)
minb2=matrix(,1000,14)
maxb2=matrix(,1000,14)
tfb2=matrix(,1000,14)
ftb2=matrix(,1000,14)
sfb2=matrix(,1000,14)
maxb2=matrix(,1000,14)

for (i in 1:1000){
  genetic2 <- matchit(DRUGP~ALLGEMP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+DEM+Duration10, data=data, method="genetic")
  genetic.data2 <- match.data(genetic2)
  
  pre.balance.g2 <- summary(genetic2)$sum.all
  post.balance.g2 <- summary(genetic2)$sum.matched
  
  genetic.imbalance2 <- imbalance(group=genetic.data2$DRUGP, data=genetic.data2, drop=c("DRUGP","BTTLDintLN","ccode","distance","weights"))
  
  pre2[i,]=(pre.balance.g2[15:2,1]-pre.balance.g2[15:2,2])/pre.balance.g2[15:2,3]
  post2[i,]=(post.balance.g2[15:2,1]-post.balance.g2[15:2,2])/post.balance.g2[15:2,3]
  meanstat2[i,]=genetic.imbalance2$tab[,1]
  L1s2[i,]=genetic.imbalance2$tab[,3]
  minb2[i,]=genetic.imbalance2$tab[,4]
  tfb2[i,]=genetic.imbalance2$tab[,5]
  ftb2[i,]=genetic.imbalance2$tab[,6]
  sfb2[i,]=genetic.imbalance2$tab[,7]
  maxb2[i,]=genetic.imbalance2$tab[,8]
  
  M6 <- lm(BTTLDintLN~ALLGEMP+DRUGP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+DEM+Duration10, data=genetic.data2)
  M6$se=abs(round(cl(genetic.data2, M6, genetic.data2$ccode), 2))[,3]
  coefsne2=c(coefsne2,M6$coefficients[3])
  coefs2=c(coefs2,round(exp(M6$coefficients),3)[3])
  tstats2=c(tstats2,M6$se[3])
  n2=c(n2,dim(genetic.data2)[1])
  print(i)
}

## Code for plot (see below)
pre2[is.infinite(pre2)] <- NA 
post2[is.infinite(post2)] <- NA

## Code for table in APPENDIX A
pre.imbalance <- imbalance(group=data$DRUGP, data=data, drop=c("DRUGP","BTTLDintLN","ccode","distance","weights"))
balance6=cbind(pre.imbalance$tab[,1],apply(meanstat2,2,mean),pre.imbalance$tab[,3],apply(L1s2,2,mean),
pre.imbalance$tab[,5],apply(tfb2,2,mean),pre.imbalance$tab[,6],apply(ftb2,2,mean),pre.imbalance$tab[,7],apply(sfb2,2,mean))
rownames(balance6) = c("Gem Production (In Zone)","Hydrocarbon Prod. (In Zone)","Hydrocarbon Prod. (Out Zone)","Mountainous Terrain","Conflict at Border","Weak Rebel Group","Equal Strength","Cold War", "Internationalized Conflict","Population","Ethnic Fractionalization","Ethnic Polarization","Democracy","Duration 10 Days or Less")
print(xtable(balance6,digit=4))

#########################################################
## Plot mean standardized imbalances (for both models) ##
#########################################################

par(mfrow=c(1,3), xpd=FALSE)
par(mar=c(5,0,4,0))
plot(x=apply(pre,2,mean,na.rm=T),
     y=2:15, pch=19, col="red",
     xlab="Model 2 - Total Combat Deaths", axes=FALSE, cex.lab=1,
     xlim=c(-1.4,1.4), ylab="")
points(x=apply(post,2,mean,na.rm=T),
       y=2:15, pch=19, col="blue")
segments(x0=apply(post,2,quantile,probs=c(.025),na.rm=T),y0=2:15,y1=2:15,x1=apply(post,2,quantile,probs=c(.975),na.rm=T),col="blue")
abline(v=0, lty="dashed")
axis(1)
axis(4, at=15:2, labels=FALSE,
     las=1,cex.axis=.5)
par(xpd=TRUE)
plot(x=c(1:14), 
     y=c(1:14),axes=FALSE,
     pch=NA, 
     text(x=7.5, y=14:1, labels = cbind("Gem Production (In Zone)","Hydrocarbon Prod. (In Zone)","Hydrocarbon Prod. (Out of Zone)","Mountainous Terrain","Conflict at Border","Weak Rebel Group","Equal Strength","Cold War", "Internationalized Conflict","Population","Ethnic Fractionalization","Ethnic Polarization","Democracy","Duration 10 Days or Less")),
     xlab="",ylab="")
legend(x=4.2,y=0,legend=c("Pre-match","Post-match"),
       pch=c(19,19),
       ncol=1,
       col=c("red","blue"),
       cex=0.9,
       bty="n")
par(xpd=FALSE)
plot(x=apply(pre2,2,mean,na.rm=T),
     y=2:15, pch=19, col="red",
     xlab="Model 6 - Combat Intensity", axes=FALSE, cex.lab=1,
     xlim=c(-1.4,1.4), ylab="")
points(x=apply(post2,2,mean,na.rm=T),
       y=2:15, pch=19, col="blue")
segments(x0=apply(post2,2,quantile,probs=c(.025),na.rm=T),y0=2:15,y1=2:15,x1=apply(post2,2,quantile,probs=c(.975),na.rm=T),col="blue")
abline(v=0, lty="dashed")
axis(1)
axis(2, at=15:2, labels=FALSE,
     las=1,cex.axis=.5)

dev.off()

###############################
###### Random Drop Tests ######
###############################

## Model 2
keep.vars <- c("ccode","BTTLDtotLN","DRUGP","ALLGEMP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")
all.data <- repdata[keep.vars]
data <- na.omit(all.data)

rcoefs=c()
rcoefsne=c()
rtstats=c()

for (i in 1:length(n)){
  rand=sample(dim(data)[1],n[i], replace=F)
  genetic.data=data[-rand,]
  
  M2 <- lm(BTTLDtotLN~ALLGEMP+DRUGP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+DEM, data=genetic.data)
  M2$se=abs(round(cl(genetic.data, M2, genetic.data$ccode), 2))[,3]
  rcoefsne=c(rcoefsne,M2$coefficients[3])
  rcoefs=c(rcoefs,round(exp(M2$coefficients),3)[3])
  rtstats=c(rtstats,M2$se[3])
  print(i)
}

## Model 6
keep.vars <- c("ccode","BTTLDintLN","DRUGP","ALLGEMP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM","Duration10")
all.data <- repdata[keep.vars]
data <- na.omit(all.data)

rcoefs2=c()
rcoefsne2=c()
rtstats2=c()

for (i in 1:length(n)){
  rand2=sample(dim(data)[1],n[i], replace=F)
  genetic.data2=data[-rand2,]
  
  M6.2 <- lm(BTTLDintLN~ALLGEMP+DRUGP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+DEM+Duration10, data=genetic.data2)
  M6.2$se=abs(round(cl(genetic.data2, M6.2, genetic.data2$ccode), 2))[,3]
  rcoefsne2=c(rcoefsne2,M6.2$coefficients[3])
  rcoefs2=c(rcoefs2,round(exp(M6.2$coefficients),3)[3])
  rtstats2=c(rtstats2,M2$se[3])
  print(i)
}

########################################################################
## Plot Matched Coefs and Tstats against Random Drop Coefs and Tstats ##
########################################################################

mean(coefs)
mean(rcoefs)

## Model 2
d=density(rcoefs)
d2=density(coefs)
dt=density(rtstats)
dt2=density(tstats)

## Model 6
d.2=density(rcoefs2)
d2.2=density(coefs2)
dt.2=density(rtstats2)
dt2.2=density(tstats2)

## Plot them
par(mfrow=c(3,2), mar=c(2,4,2,2))
plot(d2,col="red",lwd=2,xlim=c(0,1.1),ylim=c(0,5.2), xlab="", main="")
lines(d,col="blue",lwd=2)
plot(d2.2,col="red",lwd=2,xlim=c(0,1.1),ylim=c(0,5.2), xlab="", main="")
lines(d.2,col="blue",lwd=2)
par(mar=c(1,4,1,2))
plot(x=c(1:12), 
     y=c(1:12),
     pch=NA, axes=FALSE,
     text(x=6, y=c(11, 6.5, 2), cex=c(1, 1.4, 1),labels = cbind("Distribution of Exponentiated Coefficients","Model 2 - Total Combat Deaths","Distribution of T-Stats")),
     xlab="",ylab="")
plot(x=c(1:12), 
     y=c(1:12),
     pch=NA, axes=FALSE,
     text(x=6, y=c(11, 6.5, 2), cex=c(1, 1.4, 1),labels = cbind("Distribution of Exponentiated Coefficients","Model 6 - Combat Intensity","Distribution of T-Stats")),
     xlab="",ylab="")
par(mar=c(3,4,0,2))
plot(dt2,col="red",lwd=2,xlim=c(0,4),ylim=c(0,1.2),xlab="", main="")
lines(dt,col="blue",lwd=2)
plot(dt2.2,col="red",lwd=2,xlim=c(0,4),ylim=c(0,1.2),xlab="", main="")
lines(dt.2,col="blue",lwd=2)


########################################################################################
########################################################################################
################################### Section 3 ########################################## 
### Re-Evaluating the Relationship between Drug Cultivation and Intrastate Violence  ###
########################################################################################
########################################################################################


###overall####
#create function to get the fd in expected battle deaths and/or intensity
#also get quantiles

qsim=function(X,data){
#Regress data and given covariates
M2b <- lm(data[,1]~as.matrix(X), data=data) 
#Robust covariance matrix, clustered by country
rvarcov=clv(data,M2b,data$ccode)
#create design matrix
Xdes=as.matrix(cbind(1,X))

#create covariates
covs=as.matrix(apply(Xdes,2,mean))

#create drug "treamtent"
covs["DRUGP",]<-1
#create drug "control"
covs2=covs
covs2["DRUGP",]<-0

####run simulation to get fd for going from no drugs to drugs####

#create expected values by averaging out the fundamental uncertainty
#i.e., create 1000 draws of the betas and sigmas and for each of those draws,
#create 1000 values of Y for drugs and no drugs and then average over them

simdata=c()
simdata2=c()
for(i in 1:1000){
  simbetas=rmvnorm(1,mean=M2b$coefficients,sigma=rvarcov)
  simsigma=sqrt(M2b$df.residual* summary(M2b)$sigma^2 / rchisq(1, M2b$df.residual))
  simdata[i]=(mean(rnorm(n=1000,(simbetas%*%covs),(simsigma))))
  simdata2[i]=(mean(rnorm(n=1000,(simbetas%*%covs2),(simsigma))))
  }

#create vector of fd's by subtracting drugs and no drugs expected values
atesetdrug=exp(simdata)-exp(simdata2)
#take average to get fd.
atedrug=(mean(atesetdrug))
#confidence intervals
low=(quantile(atesetdrug,.025))              
high=(quantile(atesetdrug,.975))  
low90=(quantile(atesetdrug,.05))              
high90=(quantile(atesetdrug,.95)) 

print(cbind(atedrug,low,high,low90,high90))
}

#################################################
### SECTION 3.1 - Total ConFlict Battle Deaths ###
#################################################

#Use the above function to see the impact of changing the battledeath threshold
# create a sequence of battle death threshold cuts from 25 to 1000

btldt=seq(25,1000,25)
fd.est=c()
low.est=c()
high.est=c()
low90.est=c()
high90.est=c()
n.est=c()
t.est=c()

for (i in btldt){
#create dataset for the new battle death threshold (i.e. keep data with only the
#threshold or higher)
  datatest=data[exp(data$BTTLDtotLN)>i,]
#create X covariates for the new dataset
  X=datatest[,c("ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord",
                "Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")]
#get fd in battle deaths going from drugs to no drugs for each dataset  
  tester=qsim(X,datatest)
  fd.est=c(fd.est,tester[1])
  low.est=c(low.est,tester[2])
  high.est=c(high.est,tester[3])
  low90.est=c(low90.est,tester[4])
  high90.est=c(high90.est,tester[5])
  n.est=c(n.est,dim(datatest)[1])
  t.est=c(t.est,sum(datatest$DRUGP))
  print(i)
}

#create table with N and number treated at each threshold
youbet=cbind(btldt,n.est,t.est)
youbet=rbind(1,youbet)
x=matrix(,dim(youbet)[1],3)
for (i in 2:dim(youbet)[1]){
x[i,1:3]=ifelse(youbet[i,2]==youbet[i-1,2],NA,i)
}
gah=na.omit(x[,1])
uniquen=(youbet[gah,])
print(xtable(uniquen, include.rownames=F,digit=0))


#FIGURE 4

par(mar=c(5,5,2,2)) 
plot(btldt,fd.est,ylim=c(round(min(low.est),-2),round(max(high.est))),ylab="First Difference in Total Battle Deaths\n(No Drugs in Conflict Zone to Drugs in Conflict Zone)"
     ,pch=20,xlab="Battle Death Threshold",cex.lab=.8,yaxt = "n",xaxt = "n",cex.axis=.7)

lines(btldt,low.est,lty=3,col="red",lwd=2)
lines(btldt,high.est,lty=3,col="red",lwd=2)
lines(btldt,low90.est,lty=3,col="blue",lwd=2)
lines(btldt,high90.est,lty=3,col="blue",lwd=2)

axis(1, at=seq(0,1000,50), labels=seq(0,1000,50),cex.axis=.7)
axis(2, at=seq(round(min(low.est),-3),round(max(high.est),-3),1000), 
labels=seq(round(min(low.est),-3),round(max(high.est),-3),1000),cex.axis=.65)
abline(h=0,lty=1)
legend(x=20,y=-6500,legend=c("95% Confidence Interval    ","90% Confidence Interval"),
      lty=c(3,3),lwd=c(1.8,1.8), col=c("red","blue"),cex=.6)

#NOW SIMULATE the effect of simply randomly deleting as many observations as are
# deleted by the threshold experiment above by sampling 100 datasets for each threshold

#first, find number of values deleted at each threshold

index=c()
for (i in 2:length(n.est)-1){
  index[i]=(dim(data)[1]-n.est[i])
}
#get unique set of number of deleted observations
del.est=unique(index)
 
fd.rans=c()
low.rans=c()
high.rans=c()
low.rans90=c()
high.rans90=c()
n.rans=c()
t.rans=c()

#start from the 2nd observation as the first is zero
for (i in 2:length(del.est)){
 
fd.ransin=c()
low.ransin=c()
high.ransin=c()
low.ransin90=c()
high.ransin90=c()
n.ransin=c()
t.ransin=c()
#now create 100 samples of deleted observations
for (j in 1:100){

  #randomly select the number of data points deleted for each threshold
  rand=sample(dim(data)[1],del.est[i], replace=F) 
  #remove randomly selected observations from the dataset
  datatest=data[-rand,]
  X=datatest[,c("ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord",
                "Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")]
  tester=qsim(X,datatest)
  fd.ransin=c(fd.ransin,tester[1])
  low.ransin=c(low.ransin,tester[2])
  high.ransin=c(high.ransin,tester[3])
  low.ransin90=c(low.ransin90,tester[4])
  high.ransin90=c(high.ransin90,tester[5])
  n.ransin=c(n.ransin,dim(datatest)[1])
  t.ransin=c(t.ransin,sum(datatest$DRUGP))
}
  fd.rans=c(fd.rans,mean(fd.ransin))
  low.rans=c(low.rans,mean(low.ransin))
  high.rans=c(high.rans,mean(high.ransin))
  low.rans90=c(low.rans90,mean(low.ransin90))
  high.rans90=c(high.rans90,mean(high.ransin90))
  n.rans=c(n.rans,mean(n.ransin))
  t.rans=c(t.rans,round(mean(t.ransin),0))
  print(i)
}

#FIGURE 5
par(mar=c(5,5,4,2)) 
plot(1:28,fd.rans,ylim=c(round(min(low.est),-2),round(max(high.est))),
ylab="First Difference in Total Battle Deaths\n(No Drugs in Conflict Zone to Drugs in Conflict Zone)"
     ,pch=20,
xlab="Number of Observations Randomly Deleted",cex.lab=.8,yaxt = "n",xaxt = "n",cex.axis=.7)

lines(1:28,low.rans,lty=3,col="red",lwd=2)
lines(1:28,high.rans,lty=3,col="red",lwd=2)
lines(1:28,low.rans90,lty=3,col="blue",lwd=2)
lines(1:28,high.rans90,lty=3,col="blue",lwd=2)
axis(1, at=1:28, labels=del.est[-1],cex.axis=.7)
axis(2, at=seq(round(min(low.est),-3),round(max(high.est),-3),1000), 
labels=seq(round(min(low.est),-3),round(max(high.est),-3),1000),cex.axis=.65)
abline(h=0,lty=1,lwd=.1,col="black")
legend(x=1,y=-6500,legend=c("95% Confidence Interval    ","90% Confidence Interval    "),
      lty=c(3,3),lwd=c(1.8,1.8), col=c("red","blue"),cex=.6)
axis(3, at=1:28,label=uniquen[-1,1],cex.axis=.65)
mtext("Corresponding Battle Death Threshold", side=3, cex=.8,line=3)

#################################################
### SECTION 3.2 - Conflict Intensity ############
#################################################

####DURATION CONTROL####

keep.vars <- c("BTTLDintLN","ALLGEMP","DRUGP","popLN","hydrocarbonP",
               "OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar",
               "internatz","popLN","ethfrac","BTTLDtotLN","ethPOL","DEM", "ccode","durationLN")  
data <- na.omit(repdata[keep.vars])
btldt=seq(25,1000,25)
fd.estint=c()
low.estint=c()
high.estint=c()
low90.estint=c()
high90.estint=c()
n.estint=c()
t.estint=c()

for (i in btldt){
  #create dataset for the new battle death threshold (i.e. keep data with only the
  #threshold or higher)
  datatest=data[exp(data$BTTLDtotLN)>i,]
  #create X covariates for the new dataset
  X=datatest[,c("ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord",
                "Weakrebel","Equalstrenght","coldwar","internatz","durationLN","popLN","ethfrac","ethPOL","DEM")]
  #get fd in battle deaths going from drugs to no drugs for each dataset  
  tester=qsim(X,datatest)
  fd.estint=c(fd.estint,tester[1])
  low.estint=c(low.estint,tester[2])
  high.estint=c(high.estint,tester[3])
  low90.estint=c(low90.estint,tester[4])
  high90.estint=c(high90.estint,tester[5])
  n.estint=c(n.estint,dim(datatest)[1])
  t.estint=c(t.estint,sum(datatest$DRUGP))
  print(i)
}

#FIGURE 6
par(mar=c(5,5,2,2)) 
plot(btldt,fd.estint,ylim=c(-11,round(max(high.estint))),ylab="First Difference in Battle Deaths per Day\n(No Drugs in Conflict Zone to Drugs in Conflict Zone)"
     ,pch=20,xlab="Battle Death Threshold",cex.lab=.8,yaxt = "n",xaxt = "n",cex.axis=.7)

lines(btldt,low.estint,lty=3,col="red",lwd=2)
lines(btldt,high.estint,lty=3,col="red",lwd=2)
lines(btldt,low90.estint,lty=3,col="blue",lwd=2)
lines(btldt,high90.estint,lty=3,col="blue",lwd=2)

axis(1, at=seq(0,1000,50), labels=seq(0,1000,50),cex.axis=.7)
axis(2, at=seq(-11,round(max(high.estint),0),1), 
labels=seq(-11,round(max(high.estint),0),1),cex.axis=.65)
abline(h=0,lty=1)
legend(x=20,y=-9.5,legend=c("95% Confidence Interval    ","90% Confidence Interval"),
      lty=c(3,3),lwd=c(1.8,1.8), col=c("red","blue"),cex=.6)

####INTENSITY WITH NO DURATION CONTROL####

keep.vars <- c("BTTLDintLN","ALLGEMP","DRUGP","popLN","hydrocarbonP",
               "OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar",
               "internatz","popLN","ethfrac","BTTLDtotLN","ethPOL","DEM", "ccode")  
data <- na.omit(repdata[keep.vars])


fdnc.estint=c()
lownc.estint=c()
highnc.estint=c()
low90nc.estint=c()
high90nc.estint=c()
n.estintnc=c()
t.estintnc=c()

for (i in btldt){
  #create dataset for the new battle death threshold (i.e. keep data with only the
  #threshold or higher)
  datatest=data[exp(data$BTTLDtotLN)>i,]
  #create X covariates for the new dataset
  X=datatest[,c("ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord",
                "Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")]
  #get fd in battle deaths going from drugs to no drugs for each dataset  
  tester=qsim(X,datatest)
  fdnc.estint=c(fdnc.estint,tester[1])
  lownc.estint=c(lownc.estint,tester[2])
  highnc.estint=c(highnc.estint,tester[3])
  low90nc.estint=c(low90nc.estint,tester[4])
  high90nc.estint=c(high90nc.estint,tester[5])
  n.estintnc=c(n.estintnc,dim(datatest)[1])
  t.estintnc=c(t.estintnc,sum(datatest$DRUGP))
  print(i)
}

#FIGURE 7
par(mar=c(5,5,2,2)) 
plot(btldt,fdnc.estint,ylim=c(round(min(lownc.estint)),4),ylab="First Difference in Battle Deaths per Day\n(No Drugs in Conflict Zone to Drugs in Conflict Zone)"
     ,pch=20,xlab="Battle Death Threshold",cex.lab=.8,yaxt = "n",xaxt = "n",cex.axis=.7)

lines(btldt,lownc.estint,lty=3,col="red",lwd=2)
lines(btldt,highnc.estint,lty=3,col="red",lwd=2)
lines(btldt,low90nc.estint,lty=3,col="blue",lwd=2)
lines(btldt,high90nc.estint,lty=3,col="blue",lwd=2)

axis(1, at=seq(0,1000,50), labels=seq(0,1000,50),cex.axis=.7)
axis(2, at=seq(round(min(lownc.estint),0),4,1), 
labels=seq(round(min(lownc.estint),0),4,1),cex.axis=.65)
abline(h=0,lty=1)
legend(x=20,y=-9.5,legend=c("95% Confidence Interval    ","90% Confidence Interval"),
      lty=c(3,3),lwd=c(1.8,1.8), col=c("red","blue"),cex=.6)




###3.3 WEIBULL DURATION MODEL using LUJALA MODEL 7 AND 8 
keep.vars <- c("BTTLDtotLN","ongoingCONF","ALLGEMP","Duration","DRUGP","hydrocarbonP",
               "OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar",
               "internatz","popLN","ethfrac","ethPOL","terr","areaLN","DEM","ccode","durationLN")  

data <- na.omit(repdata[keep.vars])

#Create censored observations
data$ongoingCONF=1-data$ongoingCONF


  btldt=seq(25,1000,25)
  fd.est=c()
  low.est=c()
  high.est=c()
  low90.est=c()
  high90.est=c()
  n.est=c()
  t.est=c()
  
  for (i in btldt){
    #create dataset for the new battle death threshold (i.e. keep data with only the
    #threshold or higher)
    datatest=data[exp(data$BTTLDtotLN)>i,]
    
    z.out=zelig(Surv(Duration,ongoingCONF)~ALLGEMP+DRUGP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+
      Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+DEM+terr+terr:hydrocarbonP,
                model="weibull",data=datatest,robust=T)
    
    x.high=setx(z.out,DRUGP=1)
    x.low=setx(z.out,DRUGP=0)
    s.out=sim(z.out,x=x.low,x1=x.high)
    fd.est=c(fd.est,mean(s.out$qi$fd))
    low.est=c(low.est,quantile(s.out$qi$fd,.025))
    high.est=c(high.est,quantile(s.out$qi$fd,.975))
    low90.est=c(low90.est,quantile(s.out$qi$fd,.05))
    high90.est=c(high90.est,quantile(s.out$qi$fd,.95))
    print(i)
  }
  
  #FIGURE 8
  par(mar=c(5,5,2,2)) 
  plot(btldt,fd.est,ylim=c(round(min(low.est)),round(max(high.est))),ylab="First Difference in Total Conflict Duration Days\n(No Drugs in Conflict Zone to Drugs in Conflict Zone)"
       ,pch=20,xlab="Battle Death Threshold",cex.lab=.8,yaxt = "n",xaxt = "n",cex.axis=.7)
  
  lines(btldt,low.est,lty=3,col="red",lwd=2)
  lines(btldt,high.est,lty=3,col="red",lwd=2)
  lines(btldt,low90.est,lty=3,col="blue",lwd=2)
  lines(btldt,high90.est,lty=3,col="blue",lwd=2)
  axis(1, at=seq(0,1000,50), labels=seq(0,1000,50),cex.axis=.7)
  axis(2, at=seq(round(min(low.est),-2),round(max(high.est),-2),1000), 
       labels=seq(round(min(low.est),-2),round(max(high.est),-2),1000),cex.axis=.65)
  abline(h=0,lty=1)
  legend(x=20,y=-6500,legend=c("95% Confidence Interval    ","90% Confidence Interval"),
         lty=c(3,3),lwd=c(1.8,1.8), col=c("red","blue"),cex=.6)


####Weibull model sensitivity (footnote in 3.3)

#######weibull model with model 2 variables

keep.vars <- c("BTTLDtotLN","ongoingCONF","ALLGEMP","Duration","DRUGP","hydrocarbonP",
               "OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar",
               "internatz","popLN","ethfrac","ethPOL","DEM","ccode")  

data <- na.omit(repdata[keep.vars])

#Create censored observations
data$ongoingCONF=1-data$ongoingCONF


  btldt=seq(25,1000,25)
  fd.est=c()
  low.est=c()
  high.est=c()
  low90.est=c()
  high90.est=c()
  n.est=c()
  t.est=c()
  
  for (i in btldt){
    #create dataset for the new battle death threshold (i.e. keep data with only the
    #threshold or higher)
    datatest=data[exp(data$BTTLDtotLN)>i,]
    
    z.out=zelig(Surv(Duration,ongoingCONF)~ALLGEMP+DRUGP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+
      Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+DEM,
                model="weibull",data=datatest,robust=T)
    
    x.high=setx(z.out,DRUGP=1)
    x.low=setx(z.out,DRUGP=0)
    s.out=sim(z.out,x=x.low,x1=x.high)
    fd.est=c(fd.est,mean(s.out$qi$fd))
    low.est=c(low.est,quantile(s.out$qi$fd,.025))
    high.est=c(high.est,quantile(s.out$qi$fd,.975))
    low90.est=c(low90.est,quantile(s.out$qi$fd,.05))
    high90.est=c(high90.est,quantile(s.out$qi$fd,.95))
    print(i)
  }
    par(mar=c(5,5,2,2)) 
  plot(btldt,fd.est,ylim=c(round(min(low.est)),round(max(high.est))),ylab="First Difference in Total Conflict Duration Days\n(No Drugs in Conflict Zone to Drugs in Conflict Zone)"
       ,pch=20,xlab="Battle Death Threshold",cex.lab=.8,yaxt = "n",xaxt = "n",cex.axis=.7)
  
  lines(btldt,low.est,lty=3,col="red",lwd=2)
  lines(btldt,high.est,lty=3,col="red",lwd=2)
  lines(btldt,low90.est,lty=3,col="blue",lwd=2)
  lines(btldt,high90.est,lty=3,col="blue",lwd=2)
  
  axis(1, at=seq(0,1000,50), labels=seq(0,1000,50),cex.axis=.7)
  axis(2, at=seq(round(min(low.est),-2),round(max(high.est),-2),1000), 
       labels=seq(round(min(low.est),-2),round(max(high.est),-2),1000),cex.axis=.65)
  abline(h=0,lty=1)
  legend(x=20,y=-6500,legend=c("95% Confidence Interval    ","90% Confidence Interval"),
         lty=c(3,3),lwd=c(1.8,1.8), col=c("red","blue"),cex=.6)


#######weibull model with hydrocarbon*terr int, gem*terr int and drugs*terr int

keep.vars <- c("BTTLDtotLN","ongoingCONF","ALLGEMP","Duration","DRUGP","hydrocarbonP",
               "OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar",
               "internatz","popLN","ethfrac","ethPOL","DEM","ccode","terr")  

data <- na.omit(repdata[keep.vars])

#Create censored observations
data$ongoingCONF=1-data$ongoingCONF


  btldt=seq(25,1000,25)
  fd.est=c()
  low.est=c()
  high.est=c()
  low90.est=c()
  high90.est=c()
  n.est=c()
  t.est=c()
  
  for (i in btldt){
    #create dataset for the new battle death threshold (i.e. keep data with only the
    #threshold or higher)
    datatest=data[exp(data$BTTLDtotLN)>i,]
    
    z.out=zelig(Surv(Duration,ongoingCONF)~ALLGEMP+DRUGP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+
      Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+DEM+terr+DRUGP:terr+ALLGEMP:terr+hydrocarbonP:terr,
                model="weibull",data=datatest,robust=T)
    
    x.high=setx(z.out,DRUGP=1)
    x.low=setx(z.out,DRUGP=0)
    s.out=sim(z.out,x=x.low,x1=x.high)
    fd.est=c(fd.est,mean(s.out$qi$fd))
    low.est=c(low.est,quantile(s.out$qi$fd,.025))
    high.est=c(high.est,quantile(s.out$qi$fd,.975))
    low90.est=c(low90.est,quantile(s.out$qi$fd,.05))
    high90.est=c(high90.est,quantile(s.out$qi$fd,.95))
    print(i)
  }
    par(mar=c(5,5,2,2)) 
  plot(btldt,fd.est,ylim=c(round(min(low.est)),round(max(high.est))),ylab="First Difference in Total Conflict Duration Days\n(No Drugs in Conflict Zone to Drugs in Conflict Zone)"
       ,pch=20,xlab="Battle Death Threshold",cex.lab=.8,yaxt = "n",xaxt = "n",cex.axis=.7)
  
  lines(btldt,low.est,lty=3,col="red",lwd=2)
  lines(btldt,high.est,lty=3,col="red",lwd=2)
  lines(btldt,low90.est,lty=3,col="blue",lwd=2)
  lines(btldt,high90.est,lty=3,col="blue",lwd=2)
  
  axis(1, at=seq(0,1000,50), labels=seq(0,1000,50),cex.axis=.7)
  axis(2, at=seq(round(min(low.est),-2),round(max(high.est),-2),1000), 
       labels=seq(round(min(low.est),-2),round(max(high.est),-2),1000),cex.axis=.65)
  abline(h=0,lty=1)
  legend(x=20,y=-6500,legend=c("95% Confidence Interval    ","90% Confidence Interval"),
         lty=c(3,3),lwd=c(1.8,1.8), col=c("red","blue"),cex=.6)


#################################################
#################################################
################## Appendix B ###################
#################################################
#################################################


####DURATION CONTROL RANDOM####

keep.vars <- c("BTTLDintLN","ALLGEMP","DRUGP","popLN","hydrocarbonP",
               "OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar",
               "internatz","popLN","ethfrac","BTTLDtotLN","ethPOL","DEM", "ccode","durationLN")  
data <- na.omit(repdata[keep.vars])
btldt=seq(25,1000,25)
fdr.estint=c()
lowr.estint=c()
highr.estint=c()
low90r.estint=c()
high90r.estint=c()
nr.estint=c()
tr.estint=c()

for (i in 2:length(del.est)){
  #create dataset for the new battle death threshold (i.e. keep data with only the
  #threshold or higher)
  
  fd.ransin=c()
  low.ransin=c()
  high.ransin=c()
  low.ransin90=c()
  high.ransin90=c()
  n.ransin=c()
  t.ransin=c()
  #now create 100 samples of deleted observations
  for (j in 1:100){
    
    #randomly select the number of data points deleted for each threshold
    rand=sample(dim(data)[1],del.est[i], replace=F) 
    #remove randomly selected observations from the dataset
    datatest=data[-rand,]
    #create X covariates for the new dataset
    X=datatest[,c("ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord",
                  "Weakrebel","Equalstrenght","coldwar","internatz","durationLN","popLN","ethfrac","ethPOL","DEM")]
    #get fd in battle deaths going from drugs to no drugs for each dataset  
    tester=qsim(X,datatest)
    fd.ransin=c(fd.ransin,tester[1])
    low.ransin=c(low.ransin,tester[2])
    high.ransin=c(high.ransin,tester[3])
    low.ransin90=c(low.ransin90,tester[4])
    high.ransin90=c(high.ransin90,tester[5])
    n.ransin=c(n.ransin,dim(datatest)[1])
    t.ransin=c(t.ransin,sum(datatest$DRUGP))
  }
  
  fdr.estint=c(fdr.estint,mean(fd.ransin))
  lowr.estint=c(lowr.estint,mean(low.ransin))
  highr.estint=c(highr.estint,mean(high.ransin))
  low90r.estint=c(low90r.estint,mean(low.ransin90))
  high90r.estint=c(high90r.estint,mean(high.ransin90))
  nr.estint=c(nr.estint,mean(n.ransin))
  tr.estint=c(tr.estint,mean(t.ransin))
  print(i)
}

#Figure 9
par(mar=c(5,5,4,2)) 
plot(1:28,fdr.estint,ylim=c(-11,4),ylab="First Difference in Battle Deaths per Day\n(No Drugs in Conflict Zone to Drugs in Conflict Zone)"
     ,pch=20,xlab="Number of Observations Randomly Deleted",cex.lab=.8,yaxt = "n",xaxt = "n",cex.axis=.7)

lines(1:28,lowr.estint,lty=3,col="red",lwd=2)
lines(1:28,highr.estint,lty=3,col="red",lwd=2)
lines(1:28,low90r.estint,lty=3,col="blue",lwd=2)
lines(1:28,high90r.estint,lty=3,col="blue",lwd=2)

axis(1, at=1:28, labels=del.est[-1],cex.axis=.7)
axis(2, at=seq(-11,4,1), 
     labels=seq(-11,4,1),cex.axis=.65)
abline(h=0,lty=1)
legend(x=20,y=-9.5,legend=c("95% Confidence Interval    ","90% Confidence Interval"),
       lty=c(3,3),lwd=c(1.8,1.8), col=c("red","blue"),cex=.6)
axis(3, at=1:28,label=uniquen[-1,1],cex.axis=.65)
mtext("Corresponding Battle Death Threshold", side=3, cex=.8,line=3)


####NO DURATION CONTROL RANDOM####

keep.vars <- c("BTTLDintLN","ALLGEMP","DRUGP","popLN","hydrocarbonP",
               "OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar",
               "internatz","popLN","ethfrac","BTTLDtotLN","ethPOL","DEM", "ccode")  
data <- na.omit(repdata[keep.vars])
btldt=seq(25,1000,25)
fdrncr.estint=c()
lowrncr.estint=c()
highrncr.estint=c()
low90rncr.estint=c()
high90rncr.estint=c()
nrncr.estint=c()
trncr.estint=c()

for (i in 2:length(del.est)){
  #create dataset for the new battle death threshold (i.e. keep data with only the
  #threshold or higher)
  
  fd.ransin=c()
  low.ransin=c()
  high.ransin=c()
  low.ransin90=c()
  high.ransin90=c()
  n.ransin=c()
  t.ransin=c()
  #now create 100 samples of deleted observations
  for (j in 1:100){
    
    #randomly select the number of data points deleted for each threshold
    rand=sample(dim(data)[1],del.est[i], replace=F) 
    #remove randomly selected observations from the dataset
    datatest=data[-rand,]
    #create X covariates for the new dataset
    X=datatest[,c("ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord",
                  "Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")]
    #get fd in battle deaths going from drugs to no drugs for each dataset  
    tester=qsim(X,datatest)
    fd.ransin=c(fd.ransin,tester[1])
    low.ransin=c(low.ransin,tester[2])
    high.ransin=c(high.ransin,tester[3])
    low.ransin90=c(low.ransin90,tester[4])
    high.ransin90=c(high.ransin90,tester[5])
    n.ransin=c(n.ransin,dim(datatest)[1])
    t.ransin=c(t.ransin,sum(datatest$DRUGP))
  }
  
  fdrncr.estint=c(fdrncr.estint,mean(fd.ransin))
  lowrncr.estint=c(lowrncr.estint,mean(low.ransin))
  highrncr.estint=c(highrncr.estint,mean(high.ransin))
  low90rncr.estint=c(low90rncr.estint,mean(low.ransin90))
  high90rncr.estint=c(high90rncr.estint,mean(high.ransin90))
  nrncr.estint=c(nrncr.estint,mean(n.ransin))
  trncr.estint=c(trncr.estint,mean(t.ransin))
  print(i)
}

#Figure 10

par(mar=c(5,5,4,2)) 
plot(1:28,fdrncr.estint,ylim=c(-11,4),ylab="First Difference in Battle Deaths per Day\n(No Drugs in Conflict Zone to Drugs in Conflict Zone)"
     ,pch=20,xlab="Number of Observations Randomly Deleted",cex.lab=.8,yaxt = "n",xaxt = "n",cex.axis=.7)

lines(1:28,lowrncr.estint,lty=3,col="red",lwd=2)
lines(1:28,highrncr.estint,lty=3,col="red",lwd=2)
lines(1:28,low90rncr.estint,lty=3,col="blue",lwd=2)
lines(1:28,high90rncr.estint,lty=3,col="blue",lwd=2)

axis(1, at=1:28, labels=del.est[-1],cex.axis=.7)
axis(2, at=seq(-11,4,1), 
     labels=seq(-11,4,1),cex.axis=.65)
abline(h=0,lty=1)
legend(x=20,y=-9.5,legend=c("95% Confidence Interval    ","90% Confidence Interval"),
       lty=c(3,3),lwd=c(1.8,1.8), col=c("red","blue"),cex=.6)
axis(3, at=1:28,label=uniquen[-1,1],cex.axis=.65)
mtext("Corresponding Battle Death Threshold", side=3, cex=.8,line=3)



###RANDOM WEIBULL

fd.rans=c()
low.rans=c()
high.rans=c()
low.rans90=c()
high.rans90=c()
n.rans=c()
t.rans=c()

for (i in 2:length(del.est)){
  
  fd.ransin=c()
  low.ransin=c()
  high.ransin=c()
  low.ransin90=c()
  high.ransin90=c()
  n.ransin=c()
  t.ransin=c()
  #now create 100 samples of deleted observations
  for (j in 1:100){
    
    #randomly select the number of data points deleted for each threshold
    rand=sample(dim(data)[1],del.est[i], replace=F) 
    #remove randomly selected observations from the dataset
    datatest=data[-rand,]
    X=datatest[,c("ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord",
                  "Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM")]
    
    z.out=zelig(Surv(Duration,ongoingCONF)~ALLGEMP+DRUGP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+
      Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+DEM,
                model="weibull",data=datatest,robust=T)
    
    x.high=setx(z.out,DRUGP=1)
    x.low=setx(z.out,DRUGP=0)
    s.out=sim(z.out,x=x.low,x1=x.high)
    fd.ransin=c(fd.ransin,mean(s.out$qi$fd))
    low.ransin=c(low.ransin,quantile(s.out$qi$fd,.025))
    high.ransin=c(high.ransin,quantile(s.out$qi$fd,.975))
    low.ransin90=c(low.ransin90,quantile(s.out$qi$fd,.05))
    high.ransin90=c(high.ransin90,quantile(s.out$qi$fd,.95))
  }
  fd.rans=c(fd.rans,mean(fd.ransin))
  low.rans=c(low.rans,mean(low.ransin))
  high.rans=c(high.rans,mean(high.ransin))
  low.rans90=c(low.rans90,mean(low.ransin90))
  high.rans90=c(high.rans90,mean(high.ransin90))
  
  print(i)
}

#FIGURE 11

par(mar=c(5,5,2,2)) 
plot(1:28,fd.rans,ylim=c(round(min(low.est)),round(max(high.est))),
     ylab="First Difference in Total Conflict Duration Days\n(No Drugs in Conflict Zone to Drugs in Conflict Zone)"
     ,pch=20,
     xlab="Number of Observations Randomly Deleted",cex.lab=.8,yaxt = "n",xaxt = "n",cex.axis=.7)
lines(1:28,low.rans,lty=3,col="red",lwd=2)
lines(1:28,high.rans,lty=3,col="red",lwd=2)
lines(1:28,low.rans90,lty=3,col="blue",lwd=2)
lines(1:28,high.rans90,lty=3,col="blue",lwd=2)
axis(1, at=1:28, labels=del.est[-1],cex.axis=.7)
axis(2, at=seq(round(min(low.est),-2),round(max(high.est),-2),1000), 
     labels=seq(round(min(low.est),-2),round(max(high.est),-2),1000),cex.axis=.65)
abline(h=0,lty=1)
legend(x=1,y=-5000,legend=c("95% Confidence Interval           ","90% Confidence Interval        "),
       lty=c(3,3),lwd=c(1.8,1.8), col=c("red","blue"),cex=.6)

#################################################
#################################################
################## Appendix C ###################
#################################################
#################################################

## This script replicates Tables 3 & 4 from the article:
## Lujala, Paivi
## "Deadly Combat over Natural Resources: Gems, Petroleum, Drugs, and the Severity of Armed Civil Conflict"
## Journal of Conflict Resolution, Vol 53, No 1 (2009)

###################################################################################################
## Table 3
## OLS Equations for Total Combat Deaths and Daily Death Rate (Intensity) of Armed Civil Conflict, 1946-2002
## Model 1
keep.vars <- c("BTTLDtotLN","ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","ccode")
data <- na.omit(repdata[keep.vars])
M1 <- lm(BTTLDtotLN~ALLGEMP+DRUGP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord, data=data)
#get f stat by using robust clustered standard errors
M1urse=clv(data,M1,data$ccode)
fstat=waldtest(M1, vcov=M1urse)
#pull out t-scores and assign them to be the standard errors of the model to 
#faciliate the use of apsr table below
M1$se=abs(round(cl(data, M1, data$ccode), 2))[,3]
#exponentiate coefficients
M1$coefficients=round(exp(M1$coefficients),3)

## Model 2
keep.vars <- c("BTTLDtotLN", "ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM", "ccode")  
data <- na.omit(repdata[keep.vars])
M2 <- lm(BTTLDtotLN~ALLGEMP+DRUGP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+DEM, data=data)
#get f stat by using robust clustered standard errors
M2urse=clv(data,M2,data$ccode)
fstat=waldtest(M2, vcov=M2urse)
#pull out t-scores and assign them to be the standard errors of the model to 
#faciliate the use of apsr table below
M2$se=exp(abs(round(cl(data, M2, data$ccode), 2))[,3]
M2$coefficients=round(exp(M2$coefficients),3)

## Model 3
keep.vars <- c("BTTLDtotLN","ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","Weakrebel","Equalstrenght","internatz","DEM","ccode")  
data <- na.omit(repdata[keep.vars])
M3 <- lm(BTTLDtotLN~ALLGEMP+DRUGP+hydrocarbonP+OUTZhydrocarbonP+Weakrebel+Equalstrenght+internatz+DEM, data=data)
#get f stat by using robust clustered standard errors
M3urse=clv(data,M1,data$ccode)
fstat=waldtest(M3, vcov=M3urse)
 #pull out t-scores and assign them to be the standard errors of the model to 
#faciliate the use of apsr table below          
M3$se=abs(round(cl(data, M3, data$ccode), 2))[,3]
M3$coefficients=round(exp(M3$coefficients),3)

## Model 4 
keep.vars <- c("BTTLDtotLN","ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM","durationLN","ccode")  
data <- na.omit(repdata[keep.vars])
M4 <- lm(BTTLDtotLN~ALLGEMP+DRUGP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+DEM+durationLN, data=data)
#get f stat by using robust clustered standard errors
M4urse=clv(data,M4,data$ccode)
fstat=waldtest(M4, vcov=M4urse)
#pull out t-scores and assign them to be the standard errors of the model to 
#faciliate the use of apsr table below
M4$se=abs(round(cl(data, M4, data$ccode), 2))[,3]
M4$coefficients=round(exp(M4$coefficients),3)

## Model 5
keep.vars <- c("BTTLDintLN","ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM","ccode")  
data <- na.omit(repdata[keep.vars])
M5 <- lm(BTTLDintLN~ALLGEMP+DRUGP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+DEM, data=data)
#get f stat by using robust clustered standard errors
M5urse=clv(data,M5,data$ccode)
fstat=waldtest(M5, vcov=M5urse)
#pull out t-scores and assign them to be the standard errors of the model to 
#faciliate the use of apsr table below
M5$se=abs(round(cl(data, M5, data$ccode), 2))[,3]
M5$coefficients=round(exp(M5$coefficients),3)

## Model 6
keep.vars <- c("BTTLDintLN","ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM","Duration10","ccode")  
data <- na.omit(repdata[keep.vars])
M6 <- lm(BTTLDintLN~ALLGEMP+DRUGP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+DEM+Duration10, data=data)
#get f stat by using robust clustered standard errors
M6urse=clv(data,M6,data$ccode)
fstat=waldtest(M6, vcov=M1urse)
#pull out t-scores and assign them to be the standard errors of the model to 
#faciliate the use of apsr table below
M6$se=abs(round(cl(data, M6, data$ccode), 2))[,3]
          #exponentiate the coefficients
M6$coefficients=round(exp(M6$coefficients),3)

apsrtable(M1,M2,M3,M4,M5,M6,se="robust",digits=3,align="center",stars="NULL")

###################################################################################################
## Table 4
## OLS Estimates for Total Combat Deaths and Daily Death Rate of Armed Civil Conflict, 1946-2002: Interaction Effects

## Model 7 -- note: unable to replicate F-test for joint significance in this model
#note that the F signifigance test is added seperately into latex
keep.vars <- c("BTTLDtotLN","ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM","terr","IAterrXhydroP","ccode")  
data <- na.omit(repdata[keep.vars])

#unrestricted
M7 <- lm(BTTLDtotLN~ALLGEMP+DRUGP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+DEM+terr+IAterrXhydroP, data=data)
#get robust clustered varcov matrix
M7urse=clv(data,M7,data$ccode)
#get f stat by using robust clustered standard errors
fstat=waldtest(M7, vcov=M7urse)
 #pull out t-scores and assign them to be the standard errors of the model to 
#faciliate the use of apsr table below          
M7$se=abs(round(cl(data, M7, data$ccode), 2))[,3]
          #exponentiate the coefficients
M7$coefficients=round(exp(M7$coefficients),3)

#restricted
M7r <- lm(BTTLDtotLN~ALLGEMP+DRUGP+OUTZhydrocarbonP+mountain+confbord+Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+DEM+terr, data=data)   

#use the robust varcov matrix from the unrestricted model to get the fstat
waldtest(M7, M7r,vcov=M7urse)

## Model 8
keep.vars <- c("BTTLDintLN","ALLGEMP","DRUGP","hydrocarbonP","OUTZhydrocarbonP","mountain","confbord","Weakrebel","Equalstrenght","coldwar","internatz","popLN","ethfrac","ethPOL","DEM","Duration10","terr","IAterrXhydroP","ccode")  
data <- na.omit(repdata[keep.vars])
M8 <- lm(BTTLDintLN~ALLGEMP+DRUGP+hydrocarbonP+OUTZhydrocarbonP+mountain+confbord+Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+DEM+Duration10+terr+IAterrXhydroP, data=data)
          #get f stat by using robust clustered standard errors
M8urse=clv(data,M8,data$ccode)
fstat=waldtest(M8, vcov=M8urse)
          #pull out t-scores and assign them to be the standard errors of the model to 
          #faciliate the use of apsr table below
M8$se=abs(round(cl(data, M8, data$ccode), 2))[,3]
          #exponentiate the coefficients
          
M8$coefficients=round(exp(M8$coefficients),3)
          
M8r <- lm(BTTLDintLN~ALLGEMP+DRUGP+OUTZhydrocarbonP+mountain+confbord+Weakrebel+Equalstrenght+coldwar+internatz+popLN+ethfrac+ethPOL+DEM+Duration10+terr, data=data)
          
#get wald statistic
waldtest(M8, M8r,vcov=M8urse)
          
apsrtable(M7,M8,se="robust",digits=3,stars="NULL")